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The Gutzwiller approximation (GA) for Gutzwiller-projected grand canonical wave functions with 
fugacity factors is investigated in detail. Our systems in general contain inhomogeneity and local 
magnetic moments. In deriving renormalization formulae, we also derive or estimate terms of higher 
powers of intersite contractions neglected in the conventional GA. We examine several different 
constraints, i.e., local/global spin-dependent/independent particle-number conservation. Out of the 
four, the local spin-dependent constraint seems the most promising at present. An improved GA 
derived from it agrees with the variational Monte Carlo method better than the conventional GA 
does. The corrections to the conventional GA can be interpreted as two-site correlation including 
the phase difference of configurations. Furthermore, projected quasi-particle excited states are 
orthogonal to each other within the GA. Using these states, spectral weights are calculated. We 
show that asymmetry between electron addition and removal spectra can appear by taking into 
account the higher powers of the intersite contractions in the case of the d-wave superconductors 
and the Fermi sea; the addition is smaller than the removal. However, the asymmetry is quite weak 
especially near the Fermi level. In contrast, projected s-wave superconductors can have the opposite 
asymmetry (addition larger than removal) especially near the Fermi level. In addition, formulae from 
the other three constraints are also derived, which may be useful depending on purposes. 

PACS numbers: 71.10.Fd, 71.27.+a, 74.20.Rp 



I. INTRODUCTION 

This paper concerns calculation of expectation values 
using projected wave functions in inhomogeneous sys- 
tems. In order to study electronic systems with repul- 
sive on-site interactions, Gutzwiller proposed projected 
wave functions^ of the form PqI^o) with the Gutzwiller 
projection operator, 

Pg = - "^T"a) . (1) 

i 

to prohibit electron double occupancy on each site. Here, 
nia = c\^Ci^ with (cia) being the creation (annihila- 
tion) operator of site i and spin a. 

Expectation values of operators by this projected wave 
function can be evaluated by the variational Monte Carlo 
method (VMC) numerically exactly within statistical er- 
rors. However, the VMC requires lots of computational 
effort for some issues. In addition, it needs one run for 
each parameter set, whereas an analytical method can 
generate more general formulae that often provide us 
some hint to understand the system. Thus, instead of the 
VMC, an analytical approximation called the Gutzwiller 
approximation (GA) is used on occasions, i.e., 

with 1*^) = Pg|vE'^), where I'^'q) have a fixed particle 
number N. The factor is the Gutzwiller renormal- 
ization factor for the operator O. If one chooses a non- 
interacting or mean-field approximated wave function as 
I'l'^), the expectation value in the r.h.s. of Eq. ([2]) can 
be easily evaluated. The renormalization factor for the 
hopping term denoted by is smaller than unity because 



it is more difficult to hop in the presence of the strong 
on-site Coulomb repulsion between electrons; that for the 
exchange interaction denoted by g"^ is larger than unity 
because each site is more often singly occupied to avoid 
the other electrons. The GA was first introduced for the 
Hubbard model by Gutzwiller—!^, then reformulated by 
Ogawa et alA A clear description of the method has been 
given by Vollhardt^=. It was also applied to a mean-field 
theory for the t-J model by Zhang et al^ Improvements 
of the GA by taking more intersite correlations have been 
made by several authors?^"^ The GA usually produces 
qualitatively correct results although it is reported that 
there are also qualitative differences in some casesi^. 

The original formulation of the GA implicitly assumes 
that a wave function before the projection has a fixed par- 
ticle number N (in the following, we call it the "canoni- 
cal scheme"). If the particle number of a wave function 
has fluctuation (the "grand canonical scheme" ) , then the 
Gutzwiller projection reduces the particle number (see 
Appendix Such reduction of the particle number may 
arouse a question whether the GA as Eq. ([2]) is valid be- 
cause this equation seems to say that the wave functions 
before and after the projection have similar properties 
except for the double occupancy; are they similar if they 
have different particle numbers? To avoid such an un- 
clear path, Anderson and Ong^^, and Edegger et alt^ 
formulated a grand canonical GA by taking the canoni- 
cal scheme as a guide. Namely, one can force the projec- 
tion not to change the average particle number, by gluing 
to Pq a fugacity factor that compensates the particle- 
number reduction. To our knowledge, the fugacity fac- 
tor was first seen in a preliminary form in the paper by 
Yokoyama and Shiba^'^ to relate the canonical and the 
grand canonical VMC. Gebhardi^ introduced position- 
and spin-dependent fugacity factors for calculational con- 
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venience of the 1/d expansion whose d ^ oo hmit cor- 
responds to the GA. They also appear in the construc- 
tion of the gossamer superconductivity by LaughUni^. 
Then, Wang et al^^ used position-dependent but spin- 
independent fugacity factors for inhomogeneous systems. 

The fugacity factors allow us freedom to choose a re- 
lation between the particle numbers before and after 
the projection, and the renormalization depends on this 
choice. Recently, Ko et ai- pointed out that two contra- 
dictory formulae of the Gutzwiller renormalization fac- 
tors in the literature actually come from two different 
choices of the fugacity factors. That is, (i) the fugacity 
factors are determined so that the projection conserves 
the local particle density of each spin direction at each 
site, or (ii) so that the projection conserves the total par- 
ticle number for each spin direction (this is the usual 
canonical-scheme constraint). Mainly for the square lat- 
tice antiferromagnet, they used the canonical scheme, 
and introduced additional position- and spin-dependent 
fugacity factors, then calculated each renormalization 
factor as a ratio of probabilities for the physical process. 

In this paper, we examine in detail several differ- 
ent choices of fugacity factors that impose local/global 
spin-dependent /independent particle- number conserva- 
tion. We adopt the grand canonical scheme, and derive 
general formulae. Some of our formulae are different from 
those by the canonical derivation. Furthermore, correc- 
tions to the conventional GA are also estimated or de- 
rived by taking intersite correlations into account. The 
structure of the paper is as follows: Sees. [Til and lllll are de- 
voted for the case (i), and Sec. lIVI for (ii). First in Sec. [Til 
we derive renormalization of the hopping and the pairing 
amplitude, the local spin moments and the exchange in- 
teraction from the local spin-dependent constraint. We 
test the formulae of the hopping amplitude by compar- 
ing with the VMC. Physical interpretations are given for 
newly derived terms. Subsequently in Sec. IIIIl we also 
check orthogonality and excitation energies of projected 
Bogoljubov quasiparticle states, and discuss asymmetry 
between positive and negative bias spectra. Next, in 
Sec. IIVI formulae from the global spin-dependent con- 
straint are derived. The formulation there includes cases 
where the particle numbers before and after the projec- 
tion are unequal. In addition, grand canonical GAs with 
local/global spin- independent constraints are briefly dis- 
cussed in Sec. El 

In our impression, the grand canonical scheme simpli- 
fies calculation in many cases because it is free from com- 
plicated configuration counting. Furthermore, system- 
atic improvement is straightforward by including terms 
from larger clusters in the linked-cluster expansion.^ The 
formulation we use is similar to the 1/d expansion by 
Metzner and VoUhardt^^, and Gcbhardi^. The lowest- 
order theory in the uniform non-superconducting limit 
of our formulation for the case (i) is equivalent to c? — > cxd 
limit of the 1/d expansion. However, in inhomogeneous 
systems and in the presence of the second and the third 
neighbor hopping, it is not clear if 1/d is a good ex- 



pansion parameter. In addition, considering future im- 
provements of the theory, it may be difhcult to define 
terms of very high order in 1/d. Therefore, we naively 
use the linked-cluster expansion as Gutzwiller's original 
formulation^ , then expand it in a power series of intersite 
contractions and neglect high order terms. Furthermore, 
we do not adhere to making derived formulae into the 
form of Eq. 0. 

Throughout this paper, we use the following nota- 
tion: A wave function before a projection is denoted 
by l^'o) and it does not have a definite particle num- 
ber and may have some inhomogeneity in general. Then, 
the wave function after the projection is represented by 
l^*) = Pj^'o), where P is a generalized projector that 
includes fugacity factors defined later. The expectation 
values of an arbitrary operator O by these wave functions 
are denoted by 

Furthermore, 

n-ta = {nia)o, nij„ = {cl^Cja)o, Ay = (cj|Ci|)o, (4) 

(5) 

In addition, denotes the spin operator at site i. 



II. LOCAL CONSTRAINT 

The Gutzwiller projection changes electron-density 
distribution in inhomogeneous systems in general. How- 
ever, by introducing fugacity factors, one can force de- 
sired electron-density distribution. We prefer to start 
from the grand canonical GA with a local constraint for 
each spin direction, namely. 



(6) 



for any i and a. Note that this local constraint is different 
from the canonical scheme constraint that conserves the 
total particle number. However, this "local canonical" 
constraint simplifies the resultant formulae as shown in 
the following. For example, some of low order corrections 
to the GA vanish automatically. Furthermore, with this 
constraint, projected Bogoljubov quasiparticle states are 
approximately orthogonal to each other, and excitation 
energies are approximatively obtained by diagonalizing a 
renormalized Hamiltonian (shown in Sec. IIIip . 

In general, {Sf)^ and (S'f)o can be finite. Such 
cases will be discussed only in Sec. IIIDl and otherwise 
(5f)o — {Sf)o — and {cl^Cjg)o — are assumed. Fur- 
thermore, although we have d-wave superconductors in 
mind, there may be deviation from c?-wave in inhomo- 
geneous magnetic systems, and (c|^|c|^|)o (on-site pairing 
before the projection) can be non-zero. We discuss ef- 
fect of (c-|c||)o 7^ in Sec. IIII Gl and otherwise assume 
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(cJ|cJ|)o = for any i. We also do not consider triplet 
pairing of the form (c|^c^^)o and set it to zero for any i, j, 
a. The generalization to (cJg.c|g.)o ^ is straightforward. 



A. Condition for fugacity factors 

The projected wave function is defined as j^*) = P|^'o) 
with P = Yl^Pi, where 



(7) 



The local up and down particle numbers are controlled 

by X^^ , and the fugacity factors Xi^ will be determined 
later to satisfy Eq. In order to derive their explicit 
forms, let us calculate the density of cr-spin electron at 
site i, 



(8) 



In principle, by applying the Wick theorem, these ex- 
pectation values can be exactly evaluated. In practice, 
however, such calculation is quite difficult to carry out be- 
cause too many terms appear by the Wick decomposition. 
To approximate it, remember that intersite contractions, 
riijcr and , are much smaller than on-site contractions, 
Uia- An approximation to take the leading order with 
respect to the intersite contractions corresponds to the 
GA. Here, we take only on-site contractions. Then, I ^ i 
terms cancel out between the numerator and the denom- 
inator, namely, 



(9) 



5. - {Pi 



(1 - ni{){\ - n,|) + A,|n,|(l - n^) 
+A4n4(l - n^i). (10) 



Therefore, the condition to determine A^o- is given by 
Xic,(\ — nis)l'E.i = 1. By solving the simultaneous equa- 
tions for up and down spins, we obtain 



A. 



1 



(1 - - n,x) 



(11) 



The corrections to (-P^)o and {rua) can be calculated by 
taking into account intersite contractions between site i 
and other sites I ^ i. Let us calculate terms proportional 
to Inijfp. Such terms appear by the Wick decomposi- 
tion of {ni^P'^)o. We take on-site contractions for the 
sites other than i,l, and thus we only need to consider 
(nji-P;^)o. The operators in n^j — c||Cif are contracted 

with those in c|^Q| or ci^cj^ in P^. Then, the operators 
for the down spin are replaced by nj^ or 1 — n^. Namely, 
such contribution is written as 



{1 - nil) - Xi^{l - mi) + Xiinii =0 



In other words, the terms proportional to |?T-iitP vanish 
when Aio- is set as Eq. (jlip . Similarly, terms proportional 
to Afi also vanish. Therefore, with Eq. (fTT|) . we have 
(fiia) — riia + 0{njj^) + 0{Afj). Estimated corrections 
to Aicr are also of the order of njj^ or Afj . 



B. Hopping and pairing amplitude 

For the hopping term, similar calculation can be car- 
ried out. Namely, for i j, 



(4tCUc1;CjTCjxc]^ n ^I 



1 1 



Af^Aj^(l-na)(l-njx) ^ 
« (c-tCjt)o , (13) 

where we took on-site contractions except one intersite 
contraction (that is necessary) in the numerator. Then, 
the Gutzwiller renormalization factor is given by^^ 



1 - rii 



1 



1 



1 



(14) 



The next order in fact involves one more site other than 
i and j, but the second and the third order of the intersite 
contractions for such contribution vanish when Xia is set 
as Eq. pip . Therefore, the third order term involves only 
sites i and j. Namely, taking more contractions between 
i and j in Eq. P^ . 



to 
9ij] 



{l-n,i){l-nji) 



(15) 



The formula for 



is obtained by replacing as t<^4 

and Aij ^ —Aji. The rnji \niji p term in Eq. ([T5|) is from 
repulsive correlation between down-spin holes due to the 
Pauli principle: all of the four configurations in Fig. [1] 
contribute to (cJ|Cj|)o, but only (a) does to (c||Cj|). 
Then, taking into account repulsion between down-spin 
holes, (a) has less weight than the estimate by the con- 
ventional GA that neglects this correlation. 

On the other hand, the n,, i A*,- A,,; term is from su- 

(sin- 



perconducting correlation; negative for Aij — Aji 



Aji (triplet). This term 



glet), and positive for A^ 
seems related to the phase difference between the four 
configurations in Fig. [21 which appear in |^)o (before 
the projection). Our rough explanation in the case of 
riiji ~ riiji is as follows: Suppose Co, (b, Cc, Cd are co- 
efficients of the configuration (a,b,c,d) in |vI')o, and as- 
sume they are real numbers. Then, (aCd contributes to 
riiji, and — CfcCc contributes to Remember that the 
conventional GA can be derived by taking the ratio of 
the probability of configuration a^'^^ ; it implicitly assumes 
that CaCrf and — CfeCc have the same sign. Turning on the 
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superconducting correlation, configurations (a) and (b) 
as well as (c) and (d) start to correlate. Then, their con- 
tribution to the singlet order parameter before the pro- 
jection (4|c]j^ - c||c]^)o is proportional to CbCa + CcCrf- 
The magnitude of this quantity, however, is small if CaCd 
and — CbCc have the same sign. Therefore, to strengthen 
the singlet superconducting correlation, all of CaCdXbCc 
should be small. Accordingly, the weight of Fig. [TJa) 
should be smaller than the estimate by the conventional 
GA. 
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FIG. 1: Configurations contributing to (c||Cj|)o. Filled ar- 
rows represent occupied states, and open dashed arrows rep- 
resent unoccupied states. Only (a) contributes to (cj|Cj|). 



is averaged over every bond, and the statistical errors 
are negligible in the scale of this figure. For comparison 
with the GAs, is also adjusted to equalize the doping 
before the projection with that after. At small Av, the 
generalized GA agrees with the VMC very well. As Av 
increases, the deviation becomes larger. This is possibly 
because 0{Afj) term neglected in Eq. ^ may start to 
make an important contribution. 
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FIG. 3: The neatest-neighbor hopping amplitude (clo-Cja) cal- 
culated by the conventional GA [Eq. (|14p . broken line], the 
generalized GA [Eq. (fT5)l , solid line] , and the VMC [dots] for 
the projected uniform non-magnetic BCS d-wave supercon- 
ductor. 




FIG. 2: Roundabout correlation between (a) and (d) via (b) 
and (c). These configurations in I'I'o) correlate in the presence 
of superconductivity. 

We test this formula for a simplest case, i.e., the stan- 
dard uniform non-magnetic BCS d-wave superconductor. 



Uk 




(16) 



Vk 




A? 



Afc = Av (cos kx — cos ky) 



Ek ^ \lek 

^k = — 2t (cos kx + cos ky) — fi . 

The conventional GA as Eq. (|14p. the generalized GA as 
Eq. (fT5|) . and the VMC are compared in Fig. [3] for the 
nearest-neighbor hopping. Here, the generalized GA is 
done using 200x200 sites, and practically the finite-size 
effects are negligible; errors only come from neglect of the 
higher order of the intersite contractions, /i is adjusted 
to satisfy each hole concentration for each point. The 
VMC is carried out using 30x30 sites with x-antiperiodic 
y-periodic boundary condition. The hopping amplitude 



In inhomogeneous systems, there may be deviation 
from the d-wave. Since it is rather difficult to force the lo- 
cal constraint for the VMC in inhomogeneous systems, let 
us test the formula using a simpler non-c?-wave, namely, 
a uniform p-wave superconductor by redefining 



Ak = Av 



(17) 



in Eq. (|16p . The nearest-neighbor hopping amplitude in 
the x-direction is plotted in Fig. [H The generalized GA 
shows a good overall agreement with the VMC. It espe- 
cially reproduces characteristic peak at Av ~ 2t caused 
by the nijiA*jAji term in contrast to the conventional 
GA. 

The superconducting order parameters (c|^|C^|) 
calculated similarly to the hopping term, i.e., 



A* 



(18) 

The A*^ term represents the direct correlation between 
the j ti J i occupied state [Fig.ElJa)] and the empty state 
[Fig. El^d)]. The A*j|Aji|^ term contains the attractive 
correlation between holes of i | and j ii Aji is finite, i | 
and j "f tend to be simultaneously occupied or unoccu- 
pied, and it is less likely that only one of them is occupied. 
Accordingly, this effect increases weight of the configu- 
rations in Fig. [5ja) and (d), and appears as the positive 
correction in Eq. (fT8|) . The A*jnij-^n*j^ term represents 
roundabout correlation between i t and j | through i | 
and j I as depicted in Fig. O Argument similar to what 
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FIG. 4: The neatest-neighbor hopping amplitude of the 
projected "p-wave" superconductor calculated by the con- 
ventional GA [Eq. (|14|l . broken line], the generalized GA 
[Eq. CI}, solid line], and the VMC [dots]. 



1,1, and +1,-1, interchangeably, it is written as 
1 



1 r 1-1 
miTTij - ^ [(1 - "»t)(i - ^ia)(i - "iT)(i - '^ji) 

X |n»jTp(l - 2TOi)(l - 2mj)(l - - n^i) 

+ + 2m,0(l + 2mj)(l - - Uj^) 

+ |A,jf (1 - 2m,)(l + 2mj){l - n,^){l ~ n,^) 

+ |Aj,p(l + 2m,)(l - 2mj)(l - n,t)(l - 7ijx)l , (20) 



is used for the hopping amplitude (Fig. [2]) leads to the 
conclusion that the singlet correlation enhances weight 
of configurations in Fig. [5] in this case. 



' J 




FIG. 5: Roundabout correlation between (a) and (d) via (b) 
and (c) in |4'o). 



Note that Eqs. (|15|18p is mainly aimed at \i — j\ — 1. 
For next-nearest neighbors, 0{nfj) of \i — j\ = I may be 
comparable to 0{nf,j,) of — j'| = 2 and the former may 
be dominant especially in high dimensions. In general, 
as i and j separate from each other, the approximation 
by Eqs. (|15ll8p may lose accuracy. 



C. Spin moment and exchange interaction 

By definition, the local spin-z component at each site 
is not renormalized, i.e., 



(19) 



For the exchange interaction term (S^ • Sj), we take up to 
the second order of intersite contractions. Using symbols 



sfsy) 



E(4 



ith 



-RcKT<,i+A.,A;,] 



(21) 



(22) 



Here, Eq. (|^D|) seems different from what is derived as a 
ratio of probabilities for the physical process using the 
canonical scheme with the fugacity factors by Ko et ali^ 
We speculate that it possibly does not take into account 
all of the contractions above. 

To compare with the result by 1/d expansion by 
Gebhard^'', set Ay — Aji ~ and consider antiferro- 
magnets. By setting 71^^ = n-jg. in Eqs. (|20l2ip . these 
equations are reduced to 



(l-4mf)((5f5;)o-m,m, 



(1 - nj|)(l - nil) 



(1 - nii)(l - n^i) 



(23) 
(24) 



which are equivalent to the formula by the l/d 
expansioni^. However, when Aij ^ 0, renormalization of 
{SfS^) is not reduced to such a simple form, and we need 
the original formula, Eq. Note that Eqs. (|23l24p 

can be used also for non-superconducting ferromagnets. 
Namely, the local constraint leads to the conclusion that 
antiferromagnets and ferromagnets are renormalized sim- 
ilarly. This is in distinct contrast to results of the GA 
with global constraint as will be discussed in Sec. lIVFl 
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D. Systems with nonzero spin- a;?/ components 



A. Bogoljubov de Gennes (BdG) equation 



This choice of fugacity factors encounters difficulties 
when {Sf)o or {Sf)o is finite. Let us redo the derivation 
including 5± = (5±)o: 

Si = (f - ni|)(f - Uii) + Ai-f7ii-f(f - nil) 

+Kin,^il ~ n.T) + (A.T + >^^l " 1)5+5-. (26) 
The condition to determine Xi^ is 



A,. 

This is solved to give 

^ _ (f-n,t)(f-na)-5+5r 



f 



A,:, 



For a spin moment, (S*!^) = ("Sf )o, and 



(27) 

(28) 
(29) 



{Sf) - 



Hi 



n,|(f - na) + |5+P V "a(i - "^t) + 15+1 



(30) 

This renormalization factor for Sf' is larger than unity 
because it is not bound by the local constraint. Since xy 
component is renormalized differently from z component, 
approximation depends on humans' choice of z-axis. This 
asymmetry is probably related to what is discussed by Ko 
et ali^ The most reasonable choice of z-axis we think is 
making it parallel to (Si)o at each site. Then, 5+ = 
for any i. It is equivalent to formulating a GA with con- 
straints (riif -'rnii) = TLi and (S^) = (8^)0- However, such 
a GA may yield very complicated renormalization factors 
for intersite terms. One way to avoid such a complexity is 
to use spin-independent constraint as shown in Sec. IV Al 



III. LOCAL CONSTRAINT: EXCITED STATES 

The GA with the position- and spin-dependent con- 
straint discussed in the previous section has an advan- 
tage in constructing plausible excited states which are 
approximately orthogonal to each other as shown below. 

For shorthand notation, we use 

Ci = Cif , CNi^+i=c\^, (31) 

where Ni^ is the number of lattice sites. Then, the sub- 
script of this new operator runs from 1 to 2A'l, and we 
represent it by single Greek symbols as Cp. Furthermore, 
we define 



(32) 



As a preparation, let us begin with deriving a BdG 
equation by minimizing the Gutzwiller-approximated en- 
ergy following the procedure by Wang et ali^ In the fol- 
lowing, we work more on general properties of a BdG 
equation with the Gutzwiller projection, and do not use 
any Hamiltonian explicitly. However, what we have in 
mind is inhomogeneous i- J-type models. 



ija 



^ J,,S,.S, I Pg , (33) 



where the tij term with i 



j may represent local im- 
purity potentials. The zero-temperature grand potential 

il = ^^tj — M X^i CT ^io'^ "^^^ approximated by the 
GA, and represented by a function of n^^, namely, 



(34) 



We do not show the explicit form of ^Iqa because it can 
be derived straightforwardly by using the formulae in the 
previous section. In the derivation, one can choose the 
level of the approximation: If one takes only the leading 
order of the intersite contractions, formulae in the non- 
magnetic case are equivalent to those derived by Wang et 
al.i^ and Li et alJ^ If an improved Gutzwiller approxima- 
tion such as Eq. is used, a more accurate solution 
can be obtained in principle, although it may be more 
difficult to find self-consistent solutions. 

The chemical potential fi is determined to adjust the 
particle number N to satisfy N = —d^lcA/diJ.. The other 
variables are functional of ^'o and determined by mini- 
mizing JIga, 



sn 



GA 







(35) 



Assuming (^E'ol^'o) = 1^ then 

He = ('^*o|npcl*o) + (*o|«pcI'5*o) . (36) 
By combining Eqs. (|35l36p . 

dQcA = (<5*o|i^BdG|«'o) + (^'oli^Bdcl'^^'o) , (37) 

(38) 



where i?BdG — —r-^^n 



3n ' 

PC -Pi 

Then, f^GA takes an extremum when |\E'o) is an eigenstate 
of HBdG, namely, 



i^Bdcl^o) = ^Bdol^o) , 
SncA = £^BdG( (<5*o|*0> + (*o|^*0> 







(39) 
(40) 



The main differences from usual BdG Hamiltonian are 
the local renormalization factors in front of tij and 
Jij, and the effective local chemical potential terms 
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"SiMio-nio- with fiia- = -d^lQA/drii^ - ^ which come 
from nicr-dependence of the renormahzation factors. Lo- 
cal modulations of tij and Jij tend to be enhanced by 
the local renormahzation factors, and impurity potentials 
tend to be screened by the local chemical potentialsi^ 



B. Quasi- particles 



We rewrite -ffsdc in a matrix form, 



2N, 



(41) 



p,C=i 



The 2A'l X 2A^l matrix HpQ can be diagonalized using a 
unitary matrix U, namely. 



(42) 



Then, using 

7„ = ^(i7t)„pCp (S„>0), 
p 

7i = ^(C/t)„pCp (i?„<0), 
p 

the Hamiltonian is diagonalized as. 



(43) 



BdG 



(44) 



The ground state of this effective Hamiltonian is |^'o) = 
Yin 7n |0) . Suppose the ground state is well approximated 
by P\'^o)- Naively, one may assume that excited states 
are constructed by P7^|4'o). This form of excited states 
was first introduced for uniform systems by Zhang et alA 
For fugacity factors in P, we use those in the ground state 
even for the excited states. It probably correspond to as- 
suming that the quasi-particles 7„ are not very localized 
and that the change of the particle distribution is negli- 
gible. 



C. Orthogonality of the excited states 

The orthogonality of these excited states can be 
checked by expanding 7^ using Eq. (j43|) . For example, 
for En >0, E^> 0, 

{^ollnP PllM = E C/;„C/c™(cpP^4)o . (45) 
PC 

Here, we have to mind a discrepancy between creation 
and annihilation operators; 

PMa=4M^~f^^'^) , (46) 
C„ = Cia [(1 - flig) + Xiaflis-] . (47) 



Then, as the leading-order theory, we take on-site con- 
tractions except one intersite contraction. Thanks to 
Eq. (fTTj) , renormahzation factors are reduced to a simple 
form, i.e., 



Aio-(l-nig-) (1 - Wis) + Ajg-riis 



and we obtain simple results 



(c,^P^c],)o 



{CiaP'^Cjr)o 



(P'h 



(4o-4r>0 , 



(48) 

(49) 
(50) 

(51) 
(52) 



for any i,j,a,T including i = j. In general, (c|^P^c|j 
does not satisfy this relation, but we are lucky enough to 
use the off-site pairing assumption, {cl^clg)Q — 0, then 

(c^g.P^c|g.)o ~ 0, and can exclude such an exception. Our 
projected superconducting state includes the Fermi sea 
as a special case, and these relations look different at a 
sight from those derived for the Fermi sea by Fukushima 
et al.^ In fact, however, these are identical if one remem- 
bers that P contains fugacity factors. 

Using Eq. ()49ll52p . one can transform back from Cp to 
7„ to yield 



(^o|7„P^,t.l^o) 



(7"7m>0 = 6„ 



(53) 



That is, the excited states are orthogonal to each other 
within the GA. 



D. Excitation energy 

Let |0) and \n) denote the normalized ground and ex- 
cited states. 



|0)^ 



P|*n 



^7,1 1 *c 



ViMPWoi' v^(*o|7nP^7i|*o)' 
Neglecting the second order of the difference in 



(54) 



{n\Htj\n^ 



\Htj\0) 



E^(H-pcN-(0KcI0> 

pC -pC 

{n\HBdG\n) - {0\HBdG\0) ^ En . 



(55) 



Therefore, the excitation energies are approximately the 
same as eigenenergies of the effective Hamiltonian. 
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E. Density of states 

To calculate the local density of states, we need matrix 
elements, |(n|c|^|0)p and |(n|ciCT|0)p. First of all, using 
Eq. (|53p with n = m, the normalization of the excited 
states can be replaced as 

i^ohnPHl'^o) ~ (^olP'l^'o) . (56) 

Then, we expand 7„ in \n) using Eq. (|43p and use simple 
relations similar to Eqs. (j49ll52p . namely, 

J 



{c,rPclP)o 

{P')o 
{P')o 

{CjrPCiaP)o 

{P')o 




9ua (cLc«<t)o , 



gZ (CjrcL)c 



c'- ) 

lia \ JT 1(7 1 




(57) 
(58) 

(59) 
(60) 



for any a, t. These formulae are true also for i — j, 
more explicitly. 



{clPc^.P)0 1 



(P2 
^P4aP)0 

{P')o 



ka (1 



(61) 
(62) 



These relations are exact if exact Xia are used. 

Since the indices of the renormalization factor are only 
from those of the operator between two P's, we can trans- 
form back to 7„ to yield 

I(^l4l0>l' = .9*2.l(*o|7nct|vI/o)P 



9iicr \Upn\ 



(En > 0) 



|(n|cp|0)p = 5." l(*o|7nCp|vI/o)|^ 



C \Upn\ 



(En < 0) 



(63) 



(64) 



where p = (i^cr) as Eq. ([3T|l . The common renormaliza- 
tion factor gZ tells us that the positive and negative bias 
spectra are symmetric. This symmetric density of states 
is also obtained by the canonical scheme GAM^. We go 
one-step further about this point in the next subsection. 

For A{k,Lu), we need matrix elements in k- 
space, ^^'^ |("|cfc,o-|0) where Ck,a = 



Fourier transform of Eqs. (|57ll60p 



exp(ifci?i). These can be obtained by the 



F. Electron addition-removal asymmetry caused by 
higher order terms 

The conventional BCS theory tells us that the quasi- 
particle excitation spectra are symmetric between posi- 
tive and negative bias. However, local density of states of 



high-Tc superconductors measured by the STM is highly 
asymmetric and there is an argument that attributes this 
asymmetry to strong electron correlation.— Namely, elec- 
tron addition may be more difficult than electron removal 
because the injected electron may be repelled by the 
other electrons due to their strong Coulomb repulsion. 
It is controversial whether the projected quasi-particle 
states have symmetric spectra or not. The GA gives sym- 
metric spectrap2ii2^ if only quasi-particle excitation is con- 
sidered (incoherent excitations may cause asymmetry^) . 
In contrast, the spectra calculated by the VMC show 
asymmetryjiS, 

To discuss this point, here we calculate corrections to 
the results in the former sections. When these corrections 
are taken into account, the orthogonal relation, Eq. (|53p . 
may not be satisfied any more. Therefore, in the fol- 
lowing, we assume that the systems are almost uniform; 
in the uniform limit the wave number is a good quan- 
tum number due to the translational symmetry, and thus 
excited states are orthogonal. The next order correc- 
tions contain only site i and j similarly to those in the 
hopping term. We put general formulae in Appendix iBl 
and here only show a special case of n^i — — ni/2, 



Then, 



with 



A, 



{i 



j) 



^ 2 



(65) 
(66) 



Eqs. (|57ll60p are rewritten as 



{P')o 

{cnP4^P)o 



(P'h 



(^liP4,p)o 
{P'h 

{CjlPc^-fP)o 

{p')o 




(67) 
(68) 




5*°(4ic!T)o(l-«,^»,) , (69) 



5.L°(c,iC,t)o(l-HA,,) . 



(70) 



Since Aij > and > 0, the corrections are positive 
for the electron removal, and negative for the addition. 
For more careful analysis, we also need to check the nor- 
malization. Including the corrections, Eqs. (|49ll52p are 
rewritten as 

'\p,^^ « (sV*t)o (1 - a^aM , (71) 

-^jp^-^{cn4,)o{i-A,), (72) 
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(73) 



^^iig^«(c,xc.T)o(l-a.A,,). (74) 
V )o 

These corrections to the normalization are ah nega- 
tive, and they do not seem to cancel the asymmetry in 
Eqs. (|67II70[) . Therefore, these results suggest that the 
higher-order GA exhibits asymmetric spectra whose elec- 
tron addition spectra are smaller than the removal. 

This asymmetry is consistent with the VMC calcu- 
lations for excitation spectra by Chou et ali^, and for 
spectral weights by Bieri and Ivanov^ and Yang et al.^^ 
For more explicit comparison, we calculate the spectral 
weights, 



(7*|wfcp; namely, the spectra are just renormalized by 5*, 
and are symmetric as the standard BCS theory. In con- 
trast, by including the corrections to them, Z~ decreases 
and increases, which can cause the asymmetry in 
the spectra. These are consistent with the VMC 
results^^i^i. Note that Aij is finite even for Aij — 0, 
i.e., the Fermi sea also has the asymmetry, which is also 
consistent2i. Similarly to the hopping term, Eqs. (|67ll74p 
are more accurate for small |i — j|. Hence, the Fourier 
transformed results may include errors from the summa- 
tion over large |z — j|. It will be checked in the future 
studies by including higher order terms. Since this asym- 
metry appears as a deviation from the conventional GA, 
it is rather small (especially near the Fermi level), and 
does not look like what is seen in the STM experiment. 



Z+ik) ^ \{ka\cl\0)\ 



z-(fc) 



|(fc(T|c_fc^|0)P 



{lkaPclP)0 



{P'h 



{'ykaPC-ksP)a 



{P')o 



(75) 



(76) 



and show them in Fig. [5] for both the conventional and 
the generalized GA. Here, we include t' and t" in addition 
to Eq. ((T6)) for a better correspondence to the high-Tc 
superconductors. 




(0,0) 



(?r,0) 



(0,0) 



FIG. 6: (Color online) Z+{k) (blue lines) and Z~ (k) (red 
lines) of a projected d-wave superconductor by the conven- 
tional (dotted lines) and the generalized (solid lines) GA with 
t' = -0.3t, t" = 0.2t, Av = 0.15i, and 10% hole concentra- 
tion. 



In the case of the standard BCS theory, Z^ 



\Uk\ 



z- 



Then, for each fc-point below the Fermi level. 



one can find a corresponding point k' above the Fermi 
level such that Ek' — Ek, — Vk, ffc' = Ufc. Then, sum- 
mation of the contribution from these two points to the 
spectra is unity for both addition and removal spectr 
because \uk\^ + \uk'\^ = \vk\^ + \vk'\^ = \uk\^ + \vk\^ = 1. 
Accordingly the excitation spectra are symmetric. The 



results of the conventional GA are Z+ 



9 \Uk\ 



Z^ 



G. Opposite asymmetry in projected s-wave 
superconductors 

We speculate that the origin of the asymmetry may 
not be so simple as the intuition that electron addition 
may he more difficult than removal because electrons re- 
pel each other. Here, we show a counterexample against 
this simple scenario. That is to say, projected s-wave 
superconductors can have the opposite asymmetry; the 
electron addition spectra are larger than the removal. 
Such projected s-wave superconductors may be realized 
if the pairing interaction is isotropic because d-wave does 
not gain energy from diagonal Jij. Even if Jy is finite 
only for nearest neighbors, the mean-field approximation 
in very overdoped systems converges to extended s-wave 
solutions. To be more precise, this opposite asymmetry 
is related to finite on-site pairing before the projection, 
(c|^-|-c^)o 9^ and not really related to the symmetry of 
the gap. Then, even for d-wave, inhomogeneity causes 
deviation from the d-wave, and (c^|C^)o can be nonzero 
in general. Therefore, strongly disordered d-wave super- 
conductors could have similar properties. 

To take A^i ^ into account, we have to redo the 
derivation from the beginning. Then, and Xi should 
be replaced by 



|A,; 



l-Ui 

^2 



A,; 



2(l-n,)[(l-f)f 



(77) 
(78) 



In fact, this generalization makes analytical treatment 
very difficult, and in the following we take only the lead- 
ing order of the intersite contractions. Accordingly, its 
Aji — > limit corresponds to the conventional GA (not 
the generalized GA in Sec. lIIIF| . 

The most important matrix elements are 



{<iP<^P)^ 
{P')o 

{CilPCi'iP)o 



/A 



V^i / t t \ 



(79) 
(80) 
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Here, Eq. ((80)) is exact because of Pgci-^Pg = ci^Pg and 
Pgc\-^c\^ = 0. The renormalization factor in Eq. (17^ 
is obviously larger than that in Eq. (|5n|) . Hence, these 
matrix elements suggest that the asymmetry is the op- 
posite to that of d-wave. We have also calculated other 
matrix elements, and after the Fourier transform, for 
a miiform system are obtained as plotted in Fig. [71 







\ 


J 






J 


— 




r 



(0,0) (;r,0) {n,n) (0,0) 

k 



FIG. 7: (Color online) Z+(k) (blue lines) and Z' {k) (red 
lines) of a projected s-wave superconductor with t' = — 0.3t, 
t" = 0.2t Afc = 0.15t, and 10% hole concentration. 

In comparing Figs. [5] and [71 we should keep in mind 
that a better approximation is used for Fig. [51 and the 
asymmetry as in Fig. [6] does not appear in Fig. [71 Never- 
theless, the characteristic asymmetry near the Fermi sur- 
face in Fig. [71 is very strong, and may remain even in cal- 
culation with a better precision. Most likely, Eqs. (|79l80p 
mainly contribute to the asymmetry because the vicinity 
of the Fermi level changes most dramatically. 

H. Physical consideration for the asymmetries 

For the projected s-wave superconductors, the physical 
origin of the asymmetry may be understood as follows. 
We have been using terms "addition" and "removal" but 
these are in fact named from the ground state's view. If 
one takes the complex conjugate, this addition (removal) 
matrix elements can be regarded as removal from (ad- 
dition to) an excited state. Let us adopt the excited 
states' view for a while. In the s-wave BCS supercon- 
ducting state (before the projection), a Cooper pair may 
be formed more or less on-site, which is a resonance of 
the doubly occupied state and the empty state. When 
CiCT is operated to this wave function, it chooses the the 
"originally doubly occupied" state. Then, in Pcia\^)oi 
the opposite spin state, ia, is occupied with high proba- 
bility. Accordingly, it is and easy to remove id electron. 
In contrast, for Pc\^\^)o, it is impossible to add ia elec- 
tron. Finally, let us turn back to the ground state's view 
and review the arguments above. Then, the removal is 
difficult, hut the addition is easy. 

The asymmetry in the d-wave superconductors and the 
Fermi sea needs more consideration because it appears by 



higher order correlations. Fig.[8lja) shows a configuration 
in the ket |^')o contributing to (cj|Pc||P)o. The first 
term in Eq. I|68p represents direct correlation between Cj-\ 
and c]^| . The second term comes from the repulsive corre- 
lation between down holes, which reduces weight of this 
configuration. On the other hand, for {c'-^Pcj^P)Q, both 
configurations (a) and (b) in |^')o contribute. However, 
when electron density is high, (b) is dominant because 
empty sites are rare. Then, correlation between down 
holes increases the weight of (b). Since this effect ap- 
pears only at high density, the second term in Eq. (|67|) 
accompanies the factor a^. 

(%l 1*0) ('i'ol 1*0) 

j i j i j ' J 

(a)^|:| (b)J|:| 

FIG. 8: Configurations in |*o). 



IV. GLOBAL CONSTRAINT 
(QUASI-CANONICAL GA) 

In the former sections, we have used the local con- 
straint. However, if one needs the GA as an approximate 
method of the VMC, it may be preferred to require the 
usual canonical constraint, i.e., the total particle number 
constraint for each of up and down spins, 

^(n,.) = TVf - , (81) 

i 

where M^^^^" is the total number of a electrons after the 
projection. Although one takes N^^^^" — nio- in the 
usual canonical GA, the particle numbers before and af- 
ter the projection can be different in general. In fact, in 
the VMC, they are different; the particle number pro- 
jection P^aftc, is usually applied together with Pq, and 
the chemical potential of |^'o) is more like a variational 
parameter and does not control the particle number after 
the projection. In the following we use notation, 

after — Araftcr /at after — after , after ('fiO> 

n„ = 1\„ /iVL , n = + , (SZ) 

with TVl the total number of sites. Our purpose here is 
formulating a grand canonical GA that gives results of 
the canonical scheme, by imposing Eq. (|8ip . 

If the total spin moment is nonzero, it must be rea- 
sonable to choose the spin-z axis parallel to the global 
moment so that J2i{^f)o = J2i{Sf)o = 0. In that 
case, local xy components, (S'f)o, {Sf)o, may be finite 
in general. Then, similarly to the local-constraint for- 
mulation in Sec. Ill D| xy components of the local spin 
moments are renormalized differently from their z com- 
ponents. The canonical scheme condition for the total 
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spin moment restricts the spin-z renormalization factor 
to the vicinity of unity, but the other directions are free 
from it. We expect that this spin-rotational asymmet- 
ric rcnormahzation is a property from the canonical con- 
dition and should exist even in exact calculation. Fur- 
thermore, if the total spin moment is zero and local mo- 
ments point to various directions, we have no idea how 
to choose the z-axis. Here, to avoid such complexity, we 
assume {Sf)o = {Sf)o = 0, and as in the former sections 



A. Condition for fugacity factors 

To control the total particle numbers, we need a fac- 

tor in the form , namely, the fugacity factors 

Act do not have the site index. Accordingly, the pro- 
jected wave function is defined as l^") ~ -P|^'o) with 

P^aAf'^Af-^(l-n.T^a)- 

The formula for {fiia) has the same form as Eq. ([9]), 
and only Ao- is different, i.e., 



or third neighbor hopping in some systems. Using ai 
(1 — Ag.)(l — riia-) + Xariia, it is written as 



(87) 

The corrections to Xia and :^i affect only from third order 
and not relevant to the equation above. Note that 
goes to zero in the uniform limit with 7^^f'°'' = no-. 

In the uniform systems, by omitting irrelevant site in- 
dices and using R„ = n'^^*'"' / tier , we obtain 



9, 



I _ „aftcr I _ Jl^^^ _ Jl^^^ 



to 



1 - na Act 



(89) 



Exchange term with zero total spin-2 
component and A| = Aj 



{riicr) « nj, 



(83) 



Note that is still site-dependent because it contains 
local electron densities. By inserting it into Eq. ([5T|l . we 
obtain 



1 - 

Ao- — — = 



aftc 



(84) 



In inhomogeneous systems, A^ is solved numerically from 
Eq. (j84p in general. An important point of this uniform 
fugacity approach is that the local electron density is 
also renormalized as in Eq. (|83p. and (hia) ^ {nicy)Q in 
general. When A^ is solved and inserted into Eq. 
the corrections to (fiia) are of the second order of intersite 
contractions as will be explicitly shown in Sec. IIVDI 
The local spin-z component is renormalized as 



Act(1 - n,g) 



(85) 



where symbols t, i and +1,-1 are interchangeably used. 



The general formulae for the exchange interaction term 
are too lengthy to present here. Our main interest is in 
systems with zero total spin-z component and A| — A| 
which includes such as non-magnetic systems, antiferro- 
magnets, and stripes. Hence, for simplicity, we restrict 
ourselves to this case in the following except for Sec. lIVFl 
that treats ferromagnetic systems. The generalization to 
nonzero total spin-z component is straightforward but 
one has to work with more complexities. 



When Af = A 



A, Eq. 



(5f 



is reduced to 



(90) 



For the exchange interaction term (S^ • Sj), we take up 
to the second order of intersite contractions. Assuming 



that rrii is small, namely, rrii 
obtain 



0{nija) = 0(A,y), we 



(5j Sj ) 



4m, m, ■ 



(91) 



B. Hopping term 

The Gutzwiller renormalization factor of the hopping 
term is given by 



sfsy) 



S / QX QX I QV Qy\ 



(92) 



(clgCja) _ A^(l-n,g)(l-nj-g) _ 

Next order corrections to this formula involves another 
site, which may make important contribution for second 



where g^j = A (SjSj)~ . This is the result of the con- 
ventional GA. However, note that, if rrii is of the order 
of unity, terms such as nfj^rrii, Af jtrii, have about the 
same order of contribution as nf^^, Afj, and formulae 
above should be modified as derived below. 
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D. Beyond the "conventional" GA 

When TOj ~ riiCT, terms neglected in the previous 
derivation may grow, and we need to redo the derivation 
from the beginning. It is known that expectation values 
by Gutzwiller-projected states can be written in the form 



of a linked cluster expansion^ The terms we need here 
includes contribution from clusters one-site larger than 
those in the previous derivation. 

We relegate detailed derivation to Appendix [Cl and 
only show final results here. The renormalization of the 
particle densities is given by 



+ m c^P'i [(1 - nii)\nu^\'^ + ni||A,,,|" 



ait 



ni^lnail'^ + (1 - nii)\Aii\' 



(93) 



(94) 



where aia- = (1 — A)(l — njo-) + Anjo-. The formula for {fiia-) is obtained by replacing as t^^i ^.nd An —An. Then, 
the new equation to determine A is given by '^,i{hia) — N^^^" . The solution can be written as A « A'^^-' + A'-^-', where 
A*^°-' is A determined by Eq. and A*^^^ is the correction to it represented by 

;^(2) = <^ ^2^2(1 _ „^^)(1 _ n^) n,-f(l-ni|)+njx(l-ni|) 

^ i 

^srsr^ J "•iT(^-"u) + "-a(^~"'T) -(2) 



-ail 



{1 - 2n,i)\nat\^ - (1 - 2n,t)|A/,|2 
{1 ~ 2n,t)\nuf ~ (1 ~ 2n,^)\Aa\^ 



(95) 



Here, every A is replaced by A'^*'-' in the r.h.s. Using this new A, we can calculate spin terms. 



-(2) 



i^i / i=ii 



(2) (2) _ 

mil, mil = 



2Si^i 



au(|n,zTp + |A,,|2)-a/T(|n,;iP + |A,;|2) , (96) 



{StS^)-{St){S^) 



X miiTij 



^(2) 
1 + ^ 



„ „ I [o^ Oj)o ~ mirrij 1 — ■m^j — m 



(97) 



Here, A in the second order terms can be replaced by A*^°-'. Note that the contribution that involves a third site / in 
(Sf) cancels that in {S^S^) by the subtraction of (S'f ) There is no correction of this order for {SfSJ + SfSj) in 
Eq. (USD. 

Although several authors^i^i^ formulated improved canonical GAs by taking nearest-neighbor correlations similarly 
to ours, our result is different from any of them even if we neglect the second- and the third-neighbor terms. The 
origin of this discrepancy is not clear at present. 



E. Antiferromagnets 

As an explicit example, we show the formulae for the square lattice antiferromagnet. For periodic systems, we can 
restrict the summation over the site index i to only inside of the unit cell. In the presence of the antiferromagnetic 
moments, riij between second- or third-neighbor pairs may be comparable to or larger than that of the nearest neighbor 
pairs. To take into account these terms, we define nycr = x, x', x" , ^ij = ^ji = ^'i for the nearest, the second. 
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the third neighbor pairs, and assume these are real numbers. In addition, n^i = and Uj^ = ub for A-sublattice, 
^iT ~ "^A and = for B-sublattice, and then m — (jia — n-B)l2. By omitting irrelevant site indices, 



A(o) 



n^""-(l-"A)(l-riB) 
(1 - n^f*" )(n - 27iAns) ' 



ftA = (1 - A)(l - n^) + An^ , = (1 - A)(l - ns) + Xub 



(98) 



/I „aftcrA^aftcr 

n — ^nxnB 



S= (l-nA)(l-nB) + A(n-2nAnB) , S^,^)^ = -2 as + (4 + 4)^' , (99) 



- -(ai + a|)(x')' + 2 as (A')^ , S^^' = -(a^ + + 2 (A")^ , (100) 



"2nd 



4 A 



J (n - Iuaub) /„(2) „(2) „ 

■ ^ \^"n.n. + "2nd + " 



"(2) 
3rd 



(1 - nA)(l - n_B)(n - 2nAnB) 

+ [a^(l - 2nB) + aB(l - 2n^)] [ - + (A')' + (A")'] 

+ [a^(l - 2rM) + 05(1 - 1nB)\ [A^ - (x')' - {xf\ 



(101) 



(^l) = -(^l> 



A™ 



1 4Sil 4S^^^ 4S^;n , 4A(2A-l)m 



(2) 



^2 



X^+A2-(X')'-(A')'-(X")'-(A")' , (102) 



"(2) 



^2 

25 



^2(X' + A^) 



4to2 

1- ^(2A- 1) 



(103) 



where Eq. (|103p is for a nearest neighbor pair. In general, A^j ^ Aj^ may occur; in that case these equations need 
modification with a little more complexities. 



F. Ferromagnets 

Here, we show show a remarkable difference from the 
local constraint. That is, ferromagnets are renormalized 
very differently from antiferromagnets in contrast to re- 
sults by the local constraint in Sec. Ill CI For ferromag- 
netic wave functions without superconductivity, we can 
set Ay = 0, and Nf^'''' = Ej("»'t)o- Then, in the uni- 
form cases, this GA with the global constraint is equiv- 
alent to the one with the local-constraint. Therefore, by 
setting rrii = m in formulae in Sec. HIl we obtain those 
for the ferromagnets. 



(5f5| 



2 1 

m 

4 



, (1 - 2m) 



,(l + 2m)2 



(1 



1 , 1-n 



9 



(l-nt)(l-nx) 



1 



(104) 
(105) 
(106) 



In fact, the renormalization for the spin-z component 
represented by Eq. (|104p is different from the one de- 
rived by Zhang et al^ using a probability argument of 
the canonical GA, namely, (S'fS'p = g^(S'fS'J)o with 
defined by Eq. (|106p . However, we speculate that our re- 
sult is more reasonable because spin moment term is 
not renormalized; the canonical constraint prevents the 
spin-z component from growing larger, in contrast to the 
antiferromagnetic moments, which are not bound by the 
canonical constraint. It may be clearer if we take the 
limit of small m for Eq. (|104p . 



9h 

4 



(I. 



-yTl 



m2+g^((5f5|)o-m2) 



(107) 



and compare with Eq. (j9ip for antiferromagnets. 



Effect of iV^' 



Many of the formulae derived with the global constraint 
in this section contain a^, but it goes to zero in this 
ferromagnetic limit. It is consistent with no appearance 
of a„ in the local-constraint formulation. 



Projections reduce the Hilbert space. Hence, many 
wave functions may be equivalent to each other after the 
projection even if they are different before the projection. 
Here, we demonstrate it explicitly by the particle num- 
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ber projection. Let us start from uniform non-magnetic 
cases. Define two BCS states, 

\^o)^l[iuk+Vkclcl,^m , (108) 

k 

A: 

where N = '^^^n.i^. Under the particle number projec- 
tion P/v, 

Pn\%) = A^/'PwI^o) « Pwl^o) • (110) 

Namely, wave functions Pn\'^o) and PnI'^'q) are equiva- 
lent whereas \^o) and IvPq) are nonequivalent. At a sight, 
the quasi-particle excited states of these two BCS states 
look different because 'jI^ does not commute with A^^^. 
However, cj.^ and Ck-a in in fact yield the same 
state, and thus PnjI^I'^o) and PNl'k^l'^o) are equiva- 
lent, where jf}^ is a quasi-particle operator for l^'g). 

Therefore, even if the average particle number of |^'o) 
is not iVafter^ ^^^^ make that of |*o) equal to N''^'^" 
by choosing A to satisfy 



jyaftcr ^ ^ 



2~X'\vk\' 



Uk\ 



■X^\vk\ 



(111) 



then, using I'I'o), the GA can be applied with the local 
constraint {fiia) — {nicr)o- Accordingly, we can use con- 
venient properties derived in Sees. HIl and lllll 

Such a transformation to relate the global constraint 
to the local one may be possible also for inhomogeneous 
systems, but there seems to be a problem. The particle 
numbers can be controlled by fugacity factors Yiia '^a"^'^ 
as Eq. (|109p . One can choose \ia in |^'q) to satisfy 
("-io-) = (?iicr)o- Then, is replaced by A^^^c^^ as well 

as Cia is replaced by Cia (in annihilating the elec- 
tron, the fugacity factor caused by the creation should 
be canceled). Accordingly, I^Pq) is a rather strange wave 
function because the quasi-particles may not satisfy the 
fermion commutation relation. Then one may need to re- 
define a proper quasi-particle set for |\1/q). Furthermore, 
since the quasi-particle operators do not commute with 
the fugacity's operator, definition of the excited states 
depends on their order, in contrast to the uniform cases. 

We have originally speculated that the difference be- 
tween the particle numbers before and after the projec- 
tion may cause such asymmetry of the spectra as dis- 
cussed in the latter part of Sec. IIIII Let us look again at 
Eqs. (|61I62|) . Note that the electron removal matrix ele- 
ment is proportional to Uia- (density of the ia electrons), 
while the addition is to 1 — rii (density of the empty sites) . 
Nevertheless, it is compensated by the fugacity factor 
and the renormalization factors for the removal and the 
addition are the same. Then, one may speculate that 
some asymmetry may appear if we destroy this balance 
by changing the fugacity factors. However, according to 



our analysis here, the particle number difference does not 
have much effect, and it does not cause any asymmetry 
at least in the uniform systems. 



V. SPIN-INDEPENDENT CONSTRAINT 

At present, our main interest is in the GAs with spin- 
dependent constraints in the former sections because they 
seem to be more convenient to investigate antiferromag- 
nets, stripe state, impurity systems, and so on. However, 
in systems with a more complicated spin configuration, 
the GAs with spin-independent constraints may be use- 
ful. Therefore, here we work on it, but only take the 
leading order with respect to the intersite contractions. 



A. Local constraint 

A grand canonical GA with a spin-independent con- 
straint, 



+n^l) 



(112) 



was introduced by Wang et alr^ In non-magnetic cases 
(n^i ~ nil), this is identical to the GA with the spin- 
dependent constraint in Sees. |TT] and IIIII However, in 
magnetic cases, the results of these two GAs are differ- 
ent. Accordingly, the ferromagnetic homogeneous limitii 
with the spin-independent constraint is not equivalent 
to that of the canonical GA, but is reduced to the GA 
for charge-canonical spin-grand-canonical functions ex- 
plained in the next subsection. 

Since the formulae for and 5* have been already 
derived in Refs. [l6lfT7l . here we derive them in a slightly 
more general form by assuming that {Sf) and (Sf) are 
finite. By replacing Ajo- by A^ in the derivation in Sec. Ill Dl 
and using Sf = {Sf)o, we obtain 



6 



-, e»^(i-^.T)(i-«a)-|5+P, 

n.i 

rh^ 

(l-n,)K-2n,tn,x + 2|5+|2] ' 

(c^^Cj(j) ^ t / t t 

71 \ ^ij'^ ~ y 9iia9jja > 



_ ni{l - ni){l ~ n,_^f 




(Si • Sj 



9tj 



(113) 
(114) 

(115) 
(116) 

(117) 
(118) 

(119) 
(120) 
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In this case, the derivation of is rather simple because 
Si is nonzero only when it is operated to states where site 
i is singly occupied.— Note that (Si) || (Si)o is automat- 
ically satisfied, and there is no complexity appeared in 
Sec. Ill Dl However, in this formulation, projected quasi- 
particle excited states are not orthogonal for magnetic 
systems in general. 

B. Global constraint (charge— canonical 
spin— grand-canonical GA) 

It is possible that a wave function before the projection 
is an eigenstate of the total particle number, but not any 
eigenstate of the total spin. One can formulate a GA also 
for such systems. In this case, the condition to impose 
is the canonical condition for the total particle number, 
E..("-> : namely, 

A J2 "'"^^'"'''^ = N-''"" ■ (121) 

Here, A does not depend on a. The renormalization for- 
mulae for Sf = are the same as those in Sec. IIVI if 
Act is replaced by A. The generalization to Sf^ 7^ is 
straightforward. 

The results for antiferromagnets are equivalent to those 
in Sec. IIVI The limit to uniform ferromagncts with- 
out superconductivity can be taken by setting A^^ = 0, 
jyaftor _ J2i{f^ic)o ^ud dropping site indiceSj^^ 

^ (1 - ng){l - n)n ^ / n V 

" {\ — na)['n — 2n\ni) ' ' \n — 2n|n|/ 

(122) 

These are different from our quasi-canonical derivation in 
Eqs. (|104ll05p . In fact, in this spin-grand-canonical for- 
mulation, for ferromagncts is the same as that for anti- 
ferromagnets. This discrepancy is explained as follows: If 
the wave function before the projection is an eigenstate of 
the total spin-0, then the renormalization is represented 
by Eqs. (|104I105|) . If not, by Eq. (fT22)) . The Gutzwiller 
projection tends to magnify spin moments as explained 
in Appendix W\ However, in the canonical scheme, only 
xy components of the spins are allowed to be enhanced 
because of the canonical constraint. On the other hand, 
the spin-grand-canonical case is free from this constraint, 
and spins are renormalized more isotropically. 

VI. SUMMARY AND DISCUSSION 

We have derived various formulae using the grand 
canonical Gutzwiller approximation with several differ- 
ent constraints imposed by the fugacity factors for in- 
homogeneous magnetic systems. The formulation with 
the local particle number conservation yields more simple 
formulae. On the other hand, the global particle num- 
ber constraint is more convenient in comparing with the 
VMC. 



In Sees, mil we have discussed the asymmetry of the 
density of states. Although the incoherent spectra are 
not taken into account in this paper, we speculate that 
they appear at much higher energy scale than the coher- 
ence peaks. The conventional BCS theory tells us that 
the quasi-particle excitation spectra are symmetric be- 
tween positive and negative bias. In contrast, with the 
Gutzwiller projection, some asymmetry appears. One 
may think that electron addition is always more diffi- 
cult than electron removal if repulsion between electrons 
is strong. However, we doubt if such simple intuition 
works, and speculate that the asymmetry depends on 
the Hamiltonian. As a counterexample, we have shown 
that the projected s-wave superconductor may have the 
opposite asymmetry. Namely, even with the strong re- 
pulsion, the addition spectra can be larger than the re- 
moval. We could be able to consider in this way. Let us 
take two (normalized) Gutzwiller-projected wave func- 
tions, |'0) and 10). Suppose they are the ground states 
of Hamiltonians, and if^, respectively. Furthermore, 
we assume that |-0) and |<?!)) are also excited states of 
the other Hamiltonian, and iJ^, respectively. Then, 
(010^^1-0) is an electron addition matrix element to the 
ground state of 7/^. However, if one takes its complex 
conjugate, it is an electron removal matrix element to the 

ground state of H^. Note that Ki/ilcio-l^)! = K^lc-^rl''/')!- 
That is, if an electron addition matrix element for a 
Hamiltonian is small, an electron removal matrix element 
for a different Hamiltonian is also small. Therefore, the 
asymmetry is most likely determined not only by the pro- 
jection, but also by the Hamiltonian. 

Our ultimate purpose is to find good variational wave 
functions for systems with strong on-site Coulomb repul- 
sion. Once the fugacity factors are introduced, one pro- 
jected wave function can be related to a number of differ- 
ent unprojectedwetve functions, each of which is accompa- 
nied with fugacity factors of each definition. Therefore, 
one should probably choose a definition of fugacity fac- 
tors that matches their purpose. We speculate that the 
projected optimized solution is similar in any choice of 
the fugacity factors if the calculation is done accurately 
enough. 

There may be slight disagreement with the results by 
Ko et al. That is, their comparison between the VMC 
and the GA seems to say that the GA with the position- 
and spin-dependent constraint has larger errors than that 
with the global constraint. However, according to our 
estimation, the formulation with the position- and spin- 
dependent constraint has much smaller errors. 

To improve the approximation, one can use techniques 
of the series expansion method, such as the finite cluster 
method for calculating higher order terms. One can also 
use the Fade approximation for extrapolation if neces- 
sary, but maybe it is enough just to neglect small terms 
without extrapolation. Since this linked-cluster expan- 
sion can be done analytically, there is a possibility to 
minimize the energy analytically in contrast to the VMC. 
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APPENDIX A: DIFFERENCE BETWEEN 
CANONICAL AND GRAND CANONICAL 
SCHEME 



Suppose l^*^) is an eigenstate of the total particle 

7V|*^). Then, 



number N = Y.ia^^i<y '^i*^ ^l^o") = 
since [Pg,N]^Q, 



(Al) 



Namely, in the canonical scheme, Pq does not change 
the particle number. On the other hand in the grand 
canonical scheme, a wave function is not an eigen- 
state of iV, and the particle number is distributed with 
the distribution function, 

(0) ^ (^olPjvl^o) 
" (*o|*o) ' 

where P/v is an operator which projects onto terms with 
particle number N . Suppose p^^^ is sharply peaked at a 
mean value N ^ and fluctuation around it is of 0{^/l/N). 
Then, can be regarded as a wave function almost 
identical with l^*^) in the thermodynamic limit. In con- 
trast in the projected case, the average particle number 
of PgI^o) is in fact different from that of PgI^o')- When 
iV is large, an electron has more chance to meet another 
electron on a certain site. In other words, Pg excludes 
more states with large A^, and thus the peak position of 
the N distribution is shifted to a smaller value. Such 
distribution change was explicitly estimated by Edegger 
et alr^ as summarized below. The distribution function 
after the projection. 



Pn 



(^oIPgPatPgI^q) 
(^-oIPg-PgI^-o) 



(A2) 



can be related to that before by pN ~ gNp'N with 
_ ^{^o\PgPnPg\-^o) 



9n 



(^oIPtvI^'o) 



where C is a constant independent of A'^ coming from the 
normalization of the wave functions. The GA can esti- 
mate gN by the ratio of the relative sizes of the projected 
and unprojected Hilbert spaces as, 



Qn ~ C 



(A^L- A^t)!(^l- A^i)! 
(iVL - At - A^x)! 



where A^l is the number of lattice sites and A^| (A^i) is 
the number of up (down) spins. 

We here discuss renormalization of spins using 
Eq. (jA3p . Since [PcSi] = 0, if a wave function before 
the projection is an eigenstate of S^)^ and/or Sf 
with eigenvalues S[S -f 1),A/, respectively, then these 
quantities are not changed by Pg. If it is not such an 
eigenstate, Pg may change the distribution of S, M . If 
A is fixed, by changing M in Eq. (|A3[) . one can see that 
Pg excludes more states with small M. As a result, the 
most "probable" M increases. In fact, Eq. (|A3[) correctly 
reproduce the limit of fully polarized states [N^ = 0), 
which are obviously not affected by Pq. By rotating spin 
axes, and repeating this argument for the x, y directions, 
we conclude that Pg excludes more states with small S^. 
Physically, this can be explained as follows: To make 
small S^, electrons have to cancel their spin moments, 
and then they have high chance to meet each other on 
a certain site. When is large, the spin of an electron 
tend to orient the same direction as those of the other, 
then electrons prefer to stay at different sites, and not 
affected by the projection so much. 

The Gutzwiller projection makes more singly occu- 
pied sites, and local spin moments also tend to be mag- 
nified (this is probably related to increase of S^). In 
the canonical scheme, magnitude of uniform (ferromag- 
netic) moments are restricted by the canonical constraint, 
whereas non-uniform {e.g. antiferromagnetic, sinusoidal) 
moments are free from the canonical constraint and can 
be enhanced. On the other hand, in the spin-grand- 
canonical scheme, the total spin-z component does not 
have such restriction, and ferromagnetic moments can 
be also enhanced (shown in Sec. IV B|) . 



APPENDIX B: ELECTRON 
ADDITION/REMOVAL MATRIX ELEMENTS 



The general expressions of Eqs. (|67II74[) are written as 
follows using fiia = 1 - n.i„, aia = nicr(l - "i^)"^ and 



{cnPcl^P)o 



Into I „ ,1 "j^i'^jT^U 




(A3) 




(Bl) 
(B2) 

(B3) 
(B4) 
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{P')o 



{P')o 



-A* 



{P')o 



A,; 



(B5) 
(B6) 

(B7) 
(B8) 



APPENDIX C: HIGHER ORDER TERMS FOR 
THE GLOBAL CONSTRAINT 

By taking up to the second order of the intersite con- 
tractions, (P^)o is represented by 



{P')o 



"(2) 

1 + 



(Cl) 



Here, ^^J^ contains all the terms of the second order of 
intersite contractions in (P^)o- The division by Yli^i 
cancels single site contribution and simplifies the expres- 
sion. For calculating we need 



{n^^P^)o ^ A(l~n,4) 



E 



7(2) 



+ ^J^{aii [(1 - + n,t|A,,|2] 

-a,t [n^^\nul? + (1 - n,^)\Aa?] }. (C2) 



By taking the ratio between Eqs. (jCip and (|C2p . and 
neglecting fourth order terms, contribution from discon- 
nected clusters disappears. Then, we obtain renormal- 
ization of particle densities as Eq. 



To determine A*-^-* , we use the equation for 1 — {hi-\ 
fiti), namely. 



E 



"(2) 



EE-k 



{ai^ [(1 - 2n^)|n,,||2 - (1 - 271,t)|AhP] 

-fazT [(l-2n,t)|n,,j2-(l-2n,|)|A,,|2] }.(C3) 



By replacing A by A'^''^ + A^^-', zeroth order term cancels 
between the l.h.s. and the first term of the r.h.s. The later 
also contains A^^^; in fact, A*^^-' in the other terms can be 
negligible because they are multiplied to other intersite 
contractions. Regarding A*^^^ as the same order as nfj^ 
and Afj, and neglecting high order terms, Eq. (I95p is 
obtained. 
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